Integrate-and-fire-type models of the lateral superior olive

Neurons of the lateral superior olive (LSO) in the auditory brainstem play a fundamental role in binaural sound localization. Previous theoretical studies developed various types of neuronal models to study the physiological functions of the LSO. These models were usually tuned to a small set of physiological data with specific aims in mind. Therefore, it is unclear whether and how they can be related to each other, how widely applicable they are, and which model is suitable for what purposes. In this study, we address these questions for six different single-compartment integrate-and-fire (IF) type LSO models. The models are divided into two groups depending on their subthreshold responses: passive (linear) models with only the leak conductance and active (nonlinear) models with an additional low-voltage-activated potassium conductance that is prevalent among the auditory system. Each of these two groups is further subdivided into three subtypes according to the spike generation mechanism: one with simple threshold-crossing detection and voltage reset, one with threshold-crossing detection plus a current to mimic spike shapes, and one with a depolarizing exponential current for spiking. In our simulations, all six models were driven by identical synaptic inputs and calibrated with common criteria for binaural tuning. The resulting spike rates of the passive models were higher for intensive inputs and lower for temporally structured inputs than those of the active models, confirming the active function of the potassium current. Within each passive or active group, the simulated responses resembled each other, regardless of the spike generation types. These results, in combination with the analysis of computational costs, indicate that an active IF model is more suitable than a passive model for accurately reproducing temporal coding of LSO. The simulation of realistic spike shapes with an extended spiking mechanism added relatively small computational costs.


Introduction
Computational modeling plays a pivotal role in many fields of natural science, including neurobiology.Neuronal models in auditory neuroscience have revealed, for example, how our brain processes and combines signals from the two ears to form the perception of an acoustic image.One of the primary stations for such binaural information processing is the lateral superior olive (LSO), which is located in the auditory brainstem in the mammalian nervous system [1][2][3].As reviewed previously [4][5][6], a number of different models have been developed to simulate the physiological functions of LSO neurons.These models, however, were usually tuned with a small set of empirical data to meet the specific aims of each modeling study.Therefore, it is generally unclear how comparable these different models are to each other and how applicable they are to a wider range of questions, making it hard to re-use models in other contexts.In our previous study [7], we compared seven types of singlecompartment LSO neuron models and found that their simulated physiological responses to binaural stimuli largely resembled each other if they were calibrated with common criteria.Based on these results, we concluded that a user can readily choose an LSO model according to their intended purposes, and that a model with intermediate complexity may serve as a good starting point for general use [7], consistent with general advice to select a medium-sized model [8].
Among the LSO models we studied before [7], integrate-and-fire (IF) type neuron models have intermediate complexity.Their spike-generating mechanisms are simple thresholdcrossing detectors [9], but they still have a membrane potential that is compatible with conductance-based synaptic input models and that can be compared to physiological measurements.IF-type models have been used to simulate the physiological functions and dysfunctions of LSO in binaural processing [10][11][12][13].In the last two decades, a variety of nonlinear IF models have been developed and used [14][15][16].Among them, the exponential integrate-and-fire (EIF) model and its variations have been successfully used in many applications [15,17].We formerly predicted that an EIF-based LSO model would show similar responses to other IF models that we already examined in [7].In the present study, we aim to test this prediction with the same set of input stimuli and response criteria as used previously.In addition, we expand our list of IF-type LSO models to cover possible combinations of "add-ons".Namely, we start with the passive leaky IF model and make modifications to obtain six different models (Fig 1A) (see the first section of Results and Discussion for their general introduction).Our goals are to examine how similarly or differently these models respond to the same set of stimuli and how much computational cost is added by the modifications of each model.
The following anatomical and physiological properties of LSO comprise the bases of our modeling (see, e.g., [2,3,20,21], for more complete reviews).Excitatory inputs from spherical bushy cells (SBCs) in the ipsilateral anteroventral cochlear nucleus (AVCN), which are innervated by auditory nerve fibers, drive LSO neurons.In addition, inhibitory inputs from the medial nucleus of the trapezoid body (MNTB), which are innervated by globular bushy cells (GBCs) in the contralateral AVCN, suppress the LSO (Fig 1B).This binaural excitatory-inhibitory interaction is the basis of "anticoincidence detection" in LSO neurons [22][23][24][25][26], leading to their sensitivity to interaural time and level differences (ITDs/ILDs).The output spike rate of an LSO neuron varies periodically with the ITD of the envelope of amplitude-modulated (AM) sounds (Fig 1C), and changes monotonically with the ILD of pure tones (Fig 1D).The sensitivity of LSO neurons to these binaural cues is the neuronal basis of binaural sound localization [2,3,20,21].Furthermore, when stimulated by monaural AM tones, an LSO neuron varies its spike rate with different modulation frequencies, often showing a mild peak at around 100-400 Hz and gradual decay at higher frequencies (Fig 1E).Our modeling aims to replicate these monaural and binaural tuning properties.

IF-type LSO models to compare
The LSO comprises multiple types of neurons that can be anatomically and physiologically classified [21,27].Here we simulate the activity of sustained-spiking LSO neurons, as they have been most extensively investigated in vivo [1,20].The six IF models in this study all shared the fundamental equations (Eq 1 in "Integrate-and-fire models of LSO" of Materials and Methods).They were subdivided into two groups (Fig 1A ), according to their subthreshold responses: namely, passive IF and active IF models.The passive model (shown in green in all figures) was based on an RC circuit and its current-voltage (I-V) relationship was linear, while the active models (shown in blue) had an additional low-voltage-activated potassium (KLVA) conductance that led to a nonlinear I-V curve (Fig 2A).KLVA channels are expressed in many different neurons in the auditory system including LSO, and underlie precise temporal coding [28,29].Both passive and active models acted as low-pass filters, but the active models had  [7] with color modifications; original data in cats were collected by [18] for C and E and by [19] for D. https://doi.org/10.1371/journal.pone.0304832.g001higher impedances than the passive models at 40-400 Hz due to the activation of KLVA channels in the subthreshold regime around -60 mV (Fig 2B).
The group of passive models was further divided into three different models (Figs 1A and 2C), according to their spike-generation mechanisms.In the most basic type, called the passive leaky IF (PLk) model, a spike was counted when the simulated membrane potential reaches the threshold, which was followed by an immediate reset of the voltage to the resting state (Fig 2C,left).A possible modification of the model can be achieved by adding a spike-mimicking current [30,31].In the "passive IF with spike current" (PSp) model, a crossing of the threshold triggered a spike-like response with an up-and down-swing of the voltage (Fig 2C ,  middle).Another modification we considered was to replace the detection of threshold-crossing with a supralinear, exponential growth of the voltage [15] that resulted in a more realistic simulation of spike generation [32].This passive exponential IF (PEx) model presented a rapid voltage increase when sufficiently depolarized (Fig 2C, right), while its subthreshold responses near the resting potential remained indistinguishable from the other two passive models (Fig 2A -2C).
Analogously to the passive models, we examined three different models from the group of active IF models (Figs 1A and 2D).The active leaky IF (ALk) model and the "active IF with spike current" (ASp) model had a fixed threshold for spike generation and either voltage reset  and 2D) resembled physiological responses observed in previous in vitro studies of gerbil LSO [33].In our fitting of parameters, we carefully considered the following points to enable systematic comparisons among these models (see also "Response criteria and parameter selection" in Materials and Methods).First, we adopted the same set of parameters that were used for our previous study [7], in which the PLk and ASp models were already examined.The sole exception was the refractory period that was changed from 1.6 ms to 2.0 ms to better fit our target criteria (see "Effects of refractory period" in Methods).Second, if a parameter was shared by more than one model, then we tried to use the same value for all models, at least within the category of passive or active models.Because of these constraints, the only remaining free parameters were for the spike generating currents of the exponential models (PEx and AEx), which were chosen to satisfy the response criteria described below.

Common examination criteria
Since the main goal of this study is to compare various IF models of an LSO neuron, we fixed the input stage and evaluated the output of each model neuron.We used the same set of synaptic parameter values as in our previous comparative modeling study [7] and applied it to all six IF models.In short, each LSO neuron model was assumed to receive 20 excitatory and 8 inhibitory inputs, with their spike rate and degree of phase-locking depending on the modulation frequency (Fig 3A and 3B).These phase-locked inputs underlie the monaural and binaural LSO responses [26].To simulate ILD coding in LSO, we used a level-dependent, temporally uniform (non-phase-locked) input rate function (Fig 3C).
In our previous work [7], we set up common criteria for simulated monaural and binaural responses of modeled LSO neurons, based on available physiological data.We here adopted the same criteria to fit the parameters of the IF-type models.Namely, the maximum (peak) rate, the minimum (trough) rate, and their difference (depth) were examined for the time (Fig 3D ) and level-difference tuning curves (Fig 3E).As shown previously for the ASp model [7] and as confirmed for other models in the next sections, the response rates of the active models were all within the "targeted" range (bold numbers in Fig 3D -3E).In the passive models that lack the KLVA conductance, however, the peak rates were generally higher for the time-difference tuning curve and lower for the level-difference tuning curve than in the active models.We thus loosened the criteria by using the "accepted" ranges for binaural tuning curves for the passive models (non-bold numbers in Fig 3D -3E).The model fitting was done at the modulation frequency of 300 Hz (for phase-difference tuning) and at the ipsilateral level of 35 dB (for level-difference tuning).Once fitting of the parameters was completed, the preset "targeted" values for the monaural AM frequency tuning curves (Fig 3F ) were almost automatically satisfied by all models.In total, we set nine targeted ranges for monaural and binaural tuning curves and judged the physiological replicability of each model according to the number of ranges satisfied [7].

Simulated traces and binaural tuning curves
In this section, we describe the response characteristics of the six IF models driven by identical binaural inputs with time (phase) and level differences.The top two rows of Fig 4 (panels A 1a , A 2a , and A 3a ) show the simulated potential traces of the three passive IF models in response to stimuli phase-locked at 300 Hz.Regardless of their simulated spike shapes, their subthreshold responses and their spike timings generally resembled each other closely.This resemblance also held for level-difference tuning (Fig 4, panels A 1b , A 2b and A 3b ), in which the input spike trains were simulated with homogeneous Poisson process and had therefore no special temporal structures.These results confirm that the passive models produced almost identical responses when they were driven by the same inputs.Small differences in their spike rates and timings, however, still originated from the different spike generation mechanisms.In contrast to the immediate voltage reset in the PLk model, the membrane potential of the PSp model did not fully recover from the simulated hyperpolarization even when the absolute refractory period was over at 2 ms after threshold crossing (Fig 2C).In the PEx model, exact spike timing was slightly (around 0.1 ms) delayed from the threshold crossing times of the other two models (Fig 2C ), leading also to a delayed end of the absolute refractory period.Nevertheless, with the proper adjustment of the parameters, the phase and level tuning curves of the passive models were all within the pre-defined accepted ranges and were similar to each other (Fig 4B 1a -4B 3b ).
The observations in the passive models were mostly mirrored in the active models (Fig 4C 1a -4C 3b ).The three active-model subtypes presented very similar responses in phase-difference coding (Fig 4D 1a -4D 3a ) as well as in level-difference coding (Fig 4D 1b -4D 3b ).A notable discrepancy between the passive and active models, however, was the peak rate of the binaural tuning curves.The active models generated a larger number of spikes than the passive models for a favorable phase difference of -135 deg (compare Fig 4A 1a -4A 3a and 4C 1a -4C 3a ).In this case, hyperpolarization of the membrane caused by the phase-locked inhibitory inputs first led to the closing of KLVA channels, which in turn increased the membrane impedance (Fig 2A ) and amplified the voltage response to the succeeding phase-locked excitatory inputs, resulting A: Modeled input rates driven by AM tones with varied modulation frequencies (Eq 8 in Materials and Methods).B: Modeled degree of input phase-locking (measured as vector strength) driven by AM tones with varied modulation frequencies (Eq 9).C: Modeled input spike rates driven by unmodulated tones with varied intensities (Eq 10).D: Targeted response of LSO models driven by binaural AM tones with varied input time (modulation phase) differences.E: Targeted response of LSO models driven by binaural unmodulated tone with varied ILD.F: Targeted response of LSO models driven by monaural AM tones with varied modulation frequency.In D and E, the "targeted" value ranges are shown in bold and the slightly larger ranges of "accepted" values that were applied to passive models are shown in parentheses.https://doi.org/10.1371/journal.pone.0304832.g003

Fig 4. Simulated binaural tuning of LSO models. A 1a , A 2a , A 3a
: Simulated membrane potentials of the passive IF models driven by binaural AM tones with two different input phase differences.The phase differences of -135 deg and +45 deg corresponded approximately to the conditions in which the simulated spike rate becomes maximal and minimal, respectively.A 1b , A 2b , A 3b : Simulated membrane potentials of the passive IF models driven by binaural unmodulated tones with two different ILDs, for which the contralateral inhibitory inputs were either very weak (-10 dB) or relatively intense (+50 dB).B 1a -B 3b : Output rates of the passive models in response to binaural AM tones with varied input phase differences (B 1a , B 2a , B 3a ) or to binaural unmodulated tones with varied ILDs (B 1b , B 2b , B 3b ).C 1a -C 3b : Same plots as in A 1a -A 3b , but for the three active models.D 1a -D 3a : Same plots as in B 1a -B 3b , but for the three active models.All six models were driven by identical inputs to facilitate comparisons.In A 1a -A 3b and C 1a -C 3b , the horizontal solid and dotted lines mark the resting potentials and the spiking thresholds, respectively.In B 1a -B 3b and D 1a -D 3b , the bold numbers show the peak and trough rates, and the light pink and gray rectangular shadings represent the accepted and targeted ranges, respectively (see Fig 3D -3E for their definitions).In the leaky models (A 1a , A 1b , C 1a , C 1b ), vertical bars were manually added to show the timing of an output spike at each threshold crossing.https://doi.org/10.1371/journal.pone.0304832.g004 in an increased spiking probability.This KLVA-mediated enhancement of temporal signals is also depicted in the impedance plot (Fig 2B ), and was seen in modeling [34] and experimental [35] studies of auditory coincidence detectors.In response to intense excitatory inputs (e.g., ILD = -45 dB), the opening of the KLVA channels in the active models reduced the membrane impedance (Fig 2A ), thereby suppressing the depolarization of the membrane and resulting in a decreased spiking probability compared to the passive models (Fig 4A 1b -4A 3b and 4C 1b -4C 3b ).As the expression of KLVA conductance in LSO neurons is lower than in other prominent auditory coincidence detectors of MSO neurons or octopus cells [29,33,36], its simulated effects may not be as salient [37][38][39][40][41]. Nevertheless, our simulation results still suggest a considerable role of KLVA in suppressing the maximum response of level-difference tuning while amplifying the time-difference tuning.The combination of these effects enabled the active models to attain all targeted criteria for binaural tuning curves (Fig 3D -3E).

Further monaural and binaural responses
In the last section we confirmed that the simulated phase and level-difference tuning curves of all six models satisfy the physiology-based response criteria (Fig 3D -3E).Here we test the models with additional stimuli that were not used for parameter fitting.Simulated monaural AM frequency tuning curves of the passive (Fig 5A 1 -5A 3 ) and active (Fig 5B 1 -5B 3 ) models resembled each other, yielding curves within the pre-defined criteria (Fig 3F).The active models showed lower rates at modulation frequencies over 1 kHz, corresponding to a slightly steeper decrease of the membrane impedance with frequency (Fig 2B).Phase-tuning curves at a modulation frequency of 450 Hz were shallower in the passive models (Fig 5C 1 -5C 3 , orange) than in the active models (Fig 5D 1 -5D 3 , orange), because of the signal enhancement by the KLVA conductance discussed in the previous section.At 150 Hz modulation frequency, however, the difference between the passive and active models was small, with peak rates all within the range of 120-130 spikes/s (Fig 5C 1 -5D 3 , pink).This is due to the entrainment of spiking to the modulation frequency.Namely, the models elicited no more than one spike per modulation cycle, providing the active models with only marginal room for spike rate increase.The simulated level tuning curves showed a consistent compression of spike rates in the active models for all ipsilateral levels tested (Fig 5E 1 -5F 3 ).To sum, the active models were more sensitive to temporal signals than the passive models and more compressive to high intensity signals.Within each active or passive model group, the diversity in the spike generation mechanisms led only to small differences in the tuning curves.
Temporal discharge patterns of LSO neurons were characterized with their interspikeinterval (ISI) statistics.Earlier studies [42,43] found a variation of ISI histograms among LSO neurons in steady-state responses to monaural tonal stimulation.They also reported varied degrees of negative correlation of two succeeding ISIs.Zhou and Colburn [11] used an IF-type model equipped with afterhyperpolarization (AHP) conductance and showed that the observed cell-to-cell variability in ISI statistics can be simulated by adjusting the strength and time scale of AHP.In our simulations, both the passive (Fig 6A 1 -6A 3 ) and active (Fig 6B 1 -6B 3 ) IF models presented skewed ISI histograms, which were regarded as a signature of "fast chopping" neurons [42,43] and were associated with weak AHP [11].The hazard rate (i.e., spike occurrence rate as a function of time after preceding spike) showed a quick recovery (Fig 6C 1 -6C 3 and 6D 1 -6D 3 ), corresponding to the response of model IF neurons with rapid release from AHP [11].The simulated ISI histograms of the IF models with spike current (PSp and ASp) had a modest local peak at 2 ms (Fig 6A 2 and 6B 2 ), resembling some of those observed in vivo [43].These models were able to fire immediately at the end of (absolute) refractory period, because, unlike the other four models, the membrane potential was not fixed to the reset potential.The active models had an increased spike probability at an ISI of 3-4 ms (Fig 6D 1 -6D 3 ), which could be associated to the closing of KLVA channel during the refractory period.A similar hazard function was reported experimentally in cat LSO [42].The conditional mean of ISIs showed a flat distribution for both passive (Fig 6E 1 -6E 3 ) and active models (Fig 6F 1 -6F 3 ), indicating a lack of serial dependence of two neighboring intervals, and resembling the empirical data from a small subset of LSO neurons [42,43].In order to replicate a negative correlation between the intervals, an additional slow current that operates at a time scale of ~10 ms would have to be introduced in the model [11].

Computational efficiency
In addition to physiological replicability, computational efficiency is another determining factor in selecting a model.The rightmost column of Table 1 lists the relative computational cost of each model in comparison to the coincidence counting (CoC) model (Fig 1A), which is the simplest LSO model among those we tested.The CoC model has no membrane potential but simply operates on the coincident arrival counts of excitatory and inhibitory inputs, allowing extremely fast computation [7,26].Another reference for comparison, given in the last row of the model entries in Table 1, is the adjusted Wang-Colburn (AWC) model, a type of the Hodgkin-Huxley (HH) model [44] that simulates the temporal dynamics of several kinds of voltage- gated channels with differential equations.This model was the most complex and slowest in our collection of single-compartment LSO models [7].
The calculated computational costs of the six IF models compared in this study were all between the two extreme cases of the CoC and the AWC models (Table 1).The passive models were 2-6 times slower than the CoC model and 13-30 times faster than the HH-type AWC model.The active models, whose simulation involved a step-by-step calculation of the voltage gated KLVA conductance, required considerably larger computational costs than the passive models.They were about 20 times slower than the CoC model, but still about 4 times faster than the HH-type model.In comparison to the switch from a passive model to an active model, the costs for modifying the spike generation mechanisms were relatively low.Introducing exponential spike generation (PEx/AEx models) added slightly higher computation time than introducing a spike-mimicking current (PSp/ASp models).

Summary and concluding remarks
Having a selection of readily usable toolkits is beneficial in many fields of theoretical and experimental science, as they allow a user to select an approach suitable for the intended aims of the study.In this paper, we examined six different versions of IF models (Figs 1A and 2) and extended the list of available single-compartment LSO models [7] that were calibrated with common criteria.The simulated responses to monaural and binaural stimuli led us to the following conclusions.First, the active IF models were able to simulate physiological data more closely than the passive IF models, satisfying all the targeted response criteria (Figs 4 and  5; Table 1), but they required longer computation time.Second, the difference in spike generation mechanisms only marginally influenced the simulated physiological tunings and computational costs of both passive and active models (Table 1).Based on these results, in combination with our previous observations [7], we provide the following recommendations for selecting an LSO model: • IF-type models, either passive or active, are suitable options for investigating the interaction between modeled synaptic inputs and the membrane responses, including the simulation of the developing or aging auditory system in which the synaptic strengths are altered separately from the membrane properties (e.g., [10,13]).
• Unless the computation time is a critical factor to consider (e.g., large-scale simulations involving thousands of neurons), an active model is superior to a passive model regarding the physiological plausibility and replicability, especially for temporally structured signals.
• If simulating spike shapes or replicating the spike-initiation mechanism is not a goal, the simplified threshold-crossing detection in the PLk and ALk models usually suffices.
• When a user has no strong preference for model selection, the active leaky (ALk) model will be a good starting point as a "medium-sized" LSO neuron model.
• To simulate sequential dependence of spike times that acts on the scale of over a few milliseconds, an additional mechanism (e.g., a slow current) needs to be introduced.
• If only the input-output relationship of an LSO neuron is considered, without simulating detailed synaptic inputs, a user can also choose the CoC model [26], which is even simpler than IF-type models.Such applications include the simulation of auditory-brainstem responses [45] and perception of binaural signals [46].Limitations of our prior comparative modeling approach (see also Discussion in [7]) also apply to the present study.In our parameter fitting, we aimed to replicate "typical" LSO-like responses observed in various species.In physiological studies, however, activity of LSO neurons naturally shows variations even to the same stimuli.Therefore, fitting any model to a specific single neuron in a specific animal will require further tuning of the parameters.In addition, neurons in the LSO may have different membrane profiles along the tonotopic axis [36,47].A passive model may be suitable for simulating purely low-pass neurons, while an active model with an increased KLVA conductance may be required for simulating neurons that have band-pass membrane impedances.Furthermore, physiologically distinct subpopulations of neurons exist even within each region of LSO [21,27].In the present study, we simulated the activity of tonic spiking neurons, which continues to discharge in response to sustained sounds.Recent in vivo recordings, however, revealed a sharper binaural tuning of phasic (principal) neurons that responded only at the onset of acoustic stimulation, in comparison to tonic neurons [27,48].A two-compartment Hodgkin-Huxley-type model was developed to simulate the synaptic integration and spike generation of phasic-spiking LSO neurons [48].Creating an IF-version of phasic LSO model will be a subject of future investigations.Moreover, the modeled synaptic configurations (i.e., the number, strengths, time course, and fidelity of excitatory and inhibitory inputs) will also have to be adjusted, according to more detailed characterizations of synaptic connections to different LSO neurons in the future.
Fitting a complex HH-type model to a certain set of physiological data often requires substantial effort because of the large number of free parameters.IF-type models, in contrast, generally have far fewer parameters to fit and require a much shorter computation time for simulations.By separating and rearranging the ion-channel gating variables according to their time scales, an HH-type model can be reduced to a more easily tractable IF-type model, keeping its general spiking patterns intact (e.g., [49]).And with an appropriate modification, an IFtype model can replicate various types of spiking responses [17], providing further opportunities to simulate the physiological functions of LSO neurons efficiently and reliably with limited computational costs.We expect that the comparison of models in this study can be used as a guide for selecting a model for specific scientific questions and accompanying computational constraints.

Integrate-and-fire models of LSO
We used the same configurations for the simulated synaptic inputs to drive our IF-type LSO models as in our previous study [7].The names and the interrelations of the six IF-type models used in this study are shown by color in Fig 1A (passive models in green and active models in blue).All these models share the membrane equation of the form: which describes how the membrane potential V(t) varies with the sum I all of synaptic and membrane current.The excitatory and inhibitory synaptic currents (I synE and I synI ) are described in the next subsection.Both passive (PLk/PSp/PEx) and active (ALk/ASp/AEx) models have the leak current: In addition, the active models have the low-voltage-activated potassium (KLVA) current: Here, d(t) represents the fraction of open KLVA channels that varies with time according to the equation: In this equation, α d (V) = 0.5exp (+(V+50)/16) and β d (V) = 0.5exp(-(V+50)/16) are the activation and inactivation function of the KLVA channel in 1/ms, respectively [7].In the passive models, the KLVA conductance is fixed to zero.
Both passive and active exponential IF models (PEx/AEx) have an exponential spike-generating current of the form [15]: while in all other models this current is set to zero.In all models, a spike is counted when the membrane potential reaches the pre-set spike-detecting threshold V th .After crossing the threshold, the membrane potential of the leaky (PLk/ALk) and exponential (PEx/AEx) models is reset to the potential V ref and held at that level for a refractory period of T ref .
Threshold crossing in the models with spike currents (PSp/ASp) does not lead to an immediate voltage reset, but triggers a spike-mimicking current: with s being the time of threshold crossing [7,30].We used A 1 = 24 (nA), A 2 = 12 (nA), τ 1 = 0.17 (ms), and τ 2 = 0.37 (ms) for the passive IF model, and A 1 = 24 (nA), A 2 = 12 (nA), τ 1 = 0.15 (ms), and τ 2 = 0.30 (ms) for the active model.The resulting spike waveforms are shown in Fig 2C and 2D.The parameters of our IF models are summarized in Table 2. Justifications of the membrane parameters are provided in the "Response criteria and parameter selection" section.By eliminating the zero-current terms, the total current in the membrane equation can be written separately for each model as: In sum, all models have non-zero currents for the excitatory and inhibitory synaptic inputs and the membrane leak, while the other components additionally characterize the variation of IF models.

Common synaptic inputs
Here, we briefly describe the fundamental settings of the common input model used to drive the IF models.For more detailed descriptions and justifications, see [7,26].Based on available experimental data [33,50] and as in previous studies [7,26,44,45], each LSO neuron was assumed to receive synaptic inputs from 20 excitatory and 8 inhibitory afferents.Both excitatory and inhibitory input spike trains to the LSO model were simulated with an inhomogeneous Poisson process, whose intensity function was modeled with a von Mises distribution function that varies with time to mimic phase-locked spiking patterns of the input fibers [51,52].
For AM sound stimulation at a fixed level, the average intensity (arrival rate) function λ of the inhomogeneous Poisson process was assumed to depend on the modulation frequency: with f m being the modulation frequency in Hz.The degree of phase-locking, measured with vector strength (VS: [52,53]), was assumed to depend on the modulation frequency as: The specific numbers in these equations were determined based on empirical data (see [26] for further explanations).The curve shapes of these functions are shown in Fig 3A and 3B.To simulate phase-difference-coding of an LSO neuron, we varied the difference of the locking phases between the excitatory and inhibitory inputs.
To simulate ILD coding in response to non-modulated sounds, the average intensity function λ was assumed to be level-dependent: with s being the sound level in decibels (dB).This monotonically increasing intensity function had a sigmoidal shape (Fig 1C), mimicking the rate-level curves of bushy cells and MNTB fibers (see references in [7]).For ILD-tuning curves, the spiking activity of input fibers was assumed to be non-phase-locked (i.e., VS = 0), leading to homogeneous Poisson spike trains.
In the monaural AM stimulation, the inhibitory fibers were spontaneously active at a rate of 30 spikes per second with no phase-locking.Each input spike was converted into an alpha function [54] to obtain a unitary synaptic conductance: All spikes from all input fibers were summed to obtain the total synaptic conductance g tot (t), which was then multiplied with the difference between the synaptic reversal potential and the membrane potential to obtain the synaptic input current: To drive our IF models, we calculated the excitatory I synE and inhibitory I synI synaptic inputs separately.As detailed in [7], we adopted A syn = 3.5 (nS), τ syn = 0.16 (ms), and E syn = 0 (mV) for the excitatory synaptic input current I synE , and A syn = 12 (nS), τ syn = 0.32 (ms), and E syn = -75 (mV) for the inhibitory synaptic input current I synI .Membrane responses to simulated synaptic inputs are shown in Fig 2C and 2D.

Response criteria and parameter selection
In response to simulated binaural AM signals, the LSO models change their spike rates according to the time (phase) difference between the excitatory (ipsilateral) and inhibitory (contralateral) inputs (schematic graph in Fig 3D), which is termed the "phase-difference tuning".Driven by binaural non-phase-locked inputs, the models generally show a monotonic decrease of spike rate with the simulated ILD (schematic graph in Fig 3E), which is called the "level-difference tuning" here.To enable comparisons among the IF models of the present study and with other LSO models examined before, we re-used a set of criteria for these binaural tunings [7].The active models were tuned to fit the "targeted" spiking rates for binaural phase-difference tuning at a modulation frequency of 300 Hz and for binaural level-difference tuning at an ipsilateral level of +35 dB (bold numbers in Fig 3D -3E and in Table 1, bottom row).The passive models, however, did not satisfy all these criteria (see "Simulated traces and binaural tuning curves" in Results and Discussion for the explanation of underlying mechanisms), and therefore the "accepted" spike rates were instead used for the peak (max) values of the binaural tuning curves, as well as for the depth (the difference between peak and trough) of the binaural phase-tuning curves (non-bold numbers in Fig 3D -3E and in Table 1, bottom row).
The selection and fitting of the model parameters were performed in the following manner.
• For all six models, the membrane capacitance C, reset voltage V ref , and refractory period T ref were fixed to 24 pF, -60 mV, and 2 ms, respectively.
• For the remaining parameters of the PLk and ASp models, the same sets of parameters were adopted from our previous study [7].Therefore, these models had no free parameters to adjust.The resulting membrane resistance at -60 mV was comparable between both models (about 38 MO; Fig 2A and 2B).
• The other two passive models (PSp and PEx) adopted the values for the leak conductance g L and leak reversal potential E L from the PLk model.
• The other two active models (ALk and AEx) adopted the values for the leak conductance g L , leak reversal potential E L , KLVA conductance g K , and KLVA reversal potential E K from the ASp model.
• The PSp model used the same spike detection threshold as the PLk model.The spike-generating current of the PSp model was re-adjusted from the ASp model to fit a general spike shape of LSO [33] by accounting for the difference in their membrane impedances at depolarization (Fig 2A).At this step, no free parameter was left for the PSp model.
• The ALk models used the same spike detection threshold as the ASp model.At this step, no free parameter is left for the ALk model.
• For both exponential models (PEx and AEx), the spike generating conductance was fixed to 26.4 nS, which was the same as the leak conductance of the passive models.
• The remaining parameters of the exponential models were adjusted with grid search.A range from -50 to -40 mV (first with a coarse step of 0.5 ms, later with a finer step of 0.1 ms) for the threshold factor, from 1.2 mV to 4.0 mV (with a step of 0.2 mV) for the slope factor, and from -30 mV to 0 mV (with a step of 5 mV) for the spike detection threshold were used to find a parameter set to satisfy our binaural response criteria.
In our simulations, we calculated the average spike rates of the model neuron over 40 seconds for each parameter set.As in other models [7], more than one parameter combination attained the targeted ranges and produced similar results.Therefore, the values listed in Table 2 should be considered as one possible combination of the model parameters, but not as the "best fit".

Effects of the refractory period
In our initial attempts, we tried to fit the model with a refractory period T ref of 1.6 ms, which was adopted in our previous modeling studies [7,26].The resulting peak spike rates of some IF models, however, consistently exceeded the designated (targeted or accepted) range in our extensive grid search (see above).Hence we extended the duration of the refractory period to 2.0 ms to adjust the spike rates.This adjustment was performed in the typical range of refractory periods for cat LSO neurons (1.1-2.8 ms [55]).We also investigated how a change in the refractory period altered the spike rates of each IF model in response to binaural inputs.The peak (Fig 7A 1 -7A 2 ) and trough (Fig 7B 1 -7B 2 ) rates of the phase-difference tuning curves were almost insensitive to the change in refractory period.The PLk/PEx/ALk/AEx models (solid lines) showed a decrease in peak rate for T ref > 2.4 ms, with which the long refractory period starts to negatively impact the spike generation over a modulation period of 3.3 ms (corresponding to 300 Hz modulation).In the level-difference tuning curves, the increase in T ref linearly reduced the peak rate of these four models (Fig 7C 1 -7C 2 ), while the trough rates were not affected (Fig 7D 1 -7D 2 ).Therefore, the parameter T ref was used to finely tune the peak of leveldifference tuning without affecting other response criteria.In the PSp and ASp models, however, change in refractory period had almost no influences, because the recovery after each spike was dominated by the spike-mimicking current that lasted longer than the refractory period (Fig 2C and 2D).

Interspike-interval statistics
As an additional characterization of the IF models, we computed the spike interval statistics used in previous studies of LSO [42,43] (for theoretical accounts of the analyses, see, e.g., [56][57][58]).The LSO models were fed with the ipsilateral excitatory inputs that was used for binaural level-difference tuning curves; the contralateral inhibition was kept at the spontaneous level.ISI histograms were constructed with a time bin of T bin = 0.2 ms.The ISI count in each bin (termed N i , with the subscript i indicating the i-th bin) was normalized with the product of the bin size and the total count of ISIs, to have a unit of 1/s.The hazard rate function H i for each interval was calculated as: H i = N i / (T bin �S i ), with the normalization factor S i being the sum of all ISI counts for the i-th and higher bins [57].We stopped the calculation of the hazard function when the value of S i became smaller than 2% of the total count of ISIs.
For a second-order ISI statistic, we counted the number of occurrences for each pair of two neighboring intervals (ISI n , ISI n+1 ) in the simulated LSO spike trains.We used the same bin width T bin as above.For each bin of ISI n , we took an average of ISI n+1 to have a conditional mean [56].For a renewal process that has no memory of earlier spikes than the immediate prior, this conditional mean for ISI n+1 does not depend on the previous interval ISI n [56,58].We calculated the conditional means for intervals ISI n that had more than 0.5% of the total ISI counts.

Computation time and code availability
We used the forward Euler method for the numerical integration of the model equations with a time step 2 μs.We note, however, that a longer time step (up to around 10 μs) can generally be accepted for stable and accurate calculations of the IF models in the present study.
To evaluate the computational costs of each model, we calculated the total integration time of 40-second long traces and averaged it over 500 repetitions.In addition to the simulation of the six IF models, we also calculated the integration times of the coincidence counting (CoC) and the adjusted Wang-Colburn (AWC) models for further comparisons (see [7] for their definitions).To obtain relative computational costs, we normalized the integration time of each model by that of the CoC model, which was the fastest model introduced in our previous study [7].Numerical algorithms were implemented in Matlab (version 2021a) and simulations were carried out on a desktop computer (Dell Precision T1700) with 64-bit Windows 7 Professional Operating System, Intel Xeon CPU E3-1270 v3 (4 core, 3.5 GHz) and a 16 GB memory.The model code is publicly available at https://github.com/pinkbox-models.

Fig 1 .
Fig 1. LSO models and relevant physiological response properties.A: Interrelations of the LSO models.The six IF-type models examined in this study are color-coded consistently in all figures of this paper.Other single-compartment models shown in black were studied previously [7].Corresponding equations are given in Materials and Methods (Eqs 1-7).B: Schematic drawing of the LSO circuit.Excitatory inputs are shown in black; inhibitory inputs are shown in red.See text for abbreviations.C: Spike rate changes of a cat LSO neuron driven by AM tones with varied ITDs.Different lines correspond to different modulation frequencies.D: Spike rate changes of a cat LSO neuron driven by binaural unmodulated tones with varied ILDs.Different colors correspond to different ipsilateral sound levels.E: Spike rate changes of LSO neurons driven by monaural AM tones with varied modulation frequencies.Different lines show different units and different colors indicate different response types.Panels B-E are adapted from[7] with color modifications; original data in cats were collected by[18] for C and E and by[19] for D.

Fig 2 .
Fig 2. Sub-and suprathreshold model responses.A: Current-voltage relationship of the passive (green) and active (blue) LSO models.The spike generation mechanisms of all models were disabled for this panel.B: Membrane impedance of the passive (green) and active (blue) models.In A and B, the subtypes of active and passive models were not distinguished, because the subthreshold membrane responses were almost identical regardless of the spike-generation mechanisms.C: Responses of passive IF models to excitatory (upper panels) and inhibitory (lower panels) synaptic inputs, both of which were modeled with an alpha function (see Methods).D: Responses of active IF models to excitatory and inhibitory synaptic inputs.In each subpanel of C and D, simulated voltage traces for 1, 2, 4, 6, and 8 synchronized synaptic inputs are shown; solid and dotted horizontal lines indicate the resting and threshold voltage levels, respectively.The voltage upswing of the leaky models (PLk and ALk) were added manually to show spike timings, while the spike shapes of the other models resulted from their intrinsic features.https://doi.org/10.1371/journal.pone.0304832.g002

Fig 3 .
Fig 3. Settings of LSO model simulations.A:Modeled input rates driven by AM tones with varied modulation frequencies (Eq 8 in Materials and Methods).B: Modeled degree of input phase-locking (measured as vector strength) driven by AM tones with varied modulation frequencies (Eq 9).C: Modeled input spike rates driven by unmodulated tones with varied intensities (Eq 10).D: Targeted response of LSO models driven by binaural AM tones with varied input time (modulation phase) differences.E: Targeted response of LSO models driven by binaural unmodulated tone with varied ILD.F: Targeted response of LSO models driven by monaural AM tones with varied modulation frequency.In D and E, the "targeted" value ranges are shown in bold and the slightly larger ranges of "accepted" values that were applied to passive models are shown in parentheses.

Fig 5 .
Fig 5. Simulated responses of LSO models to stimuli not used for model fitting.A 1 -A 3 and B 1 -B 3 : Monaural AM-frequency-tuning curves of the passive (A 1 -A 3 ) and active (B 1 -B 3 ) IF models.Bold numbers show the peak rate and the rate at 1200 Hz.The gray shadings represent the targeted ranges (defined in Fig 3F).C 1 -C 3 and D 1 -D 3 : Binaural AM-phase-tuning curves of the passive (C 1 -C 3 ) and active (D 1 -D 3 ) IF models at three modulation frequencies of 150, 300 and 450 Hz.E 1 -E 3 and F 1 -F 3 : Binaural ILD-tuning curves of the passive (E 1 -E 3 ) and active (F 1 -F 3 ) IF models at five ipsilateral levels of 25-45 dB.https://doi.org/10.1371/journal.pone.0304832.g005

Fig 6 .
Fig 6.Simulated interspike-interval (ISI) statistics of LSO models to monaural tonal stimuli.A 1 -A 3 and B 1 -B 3 : ISI histograms of the passive (A 1 -A 3 ) and active (B 1 -B 3 ) IF models at two intensities (25 and 35 dB SPL).C 1 -C 3 and D 1 -D 3 : Recovery of hazard rate functions calculated from the corresponding ISI histograms.E 1 -E 3 and F 1 -F 3 : Conditional mean of ISI with respect to the preceding ISI (see Materials and Methods for the definitions of these measures).https://doi.org/10.1371/journal.pone.0304832.g006

Fig 7 .
Fig 7. Effects of refractory period on binaural tuning curves.A 1 -A 2 and B 1 -B 2 : Effect of the refractory period of the passive (A 1 and B 1 ) and active (A 2 and B 2 ) IF models on the peak (max: A 1 and A 2 ) and trough (min: B 1 and B 2 ) rates of binaural AM-phase-tuning curve at 300 Hz.C 1 -C 2 and D 1 -D 2 : Effect of the refractory period of the passive (C 1 and D 1 ) and active (C 2 and D 2 ) IF models on the peak (max: C 1 and C 2 ) and trough (min: D 1 and D 2 ) rates of binaural ILD-tuning curve at the ipsilateral level of 35 dB.The vertical dotted line shows the refractory period of 2.0 ms used for our series of simulations.The targeted and accepted ranges are indicated by the gray and light pink rectangular shadings, respectively.https://doi.org/10.1371/journal.pone.0304832.g007

Table 1 .
Summary of model output measures.Simulated spike rates that were within the targeted range are shown in bold, while those out of the targeted range but still within the accepted ranges are shown with regular fonts.Relative computation time is normalized to the average computation time of the coincidence counting model.The numbers in brackets with each model name indicate how many of the nine targeted ranges were attained.See Materials and Methods for the definition of each model and output measures.

Table 2 . Parameters of the IF models.
See text for the description of each parameter and Fig 1A for the abbreviations of the model names. https://doi.org/10.1371/journal.pone.0304832.t002